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ABSTllACT 

Three forms of the airplane spin equations of motion, derived by 
Buehler in Reference [1], fo rm the basis for the development of a 
computer program designed to seek dynamically stable equilibrium 
solutions of a spinning aircraft. The program incorporates two 
solution techniques: one based upon Euler integration, the other, a 

version of minimization by gradient search. Secondary programs are 
developed to (1) generate power-off glide parameters for use in the 
validation of the equations of motion, and (2) evaluate equation residuals 
obtained from a grid of initial conditions over the potential solution 
space. F-111 and F-4 aerodynamic force and moment models were 
utilized to evaluate the solution methods and equations of motion. The 
numerical results indicate that the F-111 and F-4 data are not repre- 
sentative of the actual aircraft and, therefore, it is highly unlikely 
that dynamically stable equilibrium solutions can be achieved from 
these models. The utility of the two solution methods is evaluated and 
the numerical results are analyzed in order to gain insight into the 
optimal application of the three forms of the equations of motion. The 
paper concludes with a discussion concerning the qualitative validation 
of the equations of motion. 
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TABLE OF SYMBOLS 



The definitions of the symbols used in this paper are as follows 



c 

g 


Center of gravity 


c 


Mean aerodynamic chord, feet 


'g 


2 

Acceleration of gravity, ft/sec 


I , I , I 

X y z 


2 

Airplane principal moments of inertia, slug~ft 


M 


Airplane mass, slugs 


N 


Load factor; number of ’’g’s” 


R 


Spin Radius, feet 


S 


Criterion function value; equal to the sum of the 
absolute values of the equation residuals, SUM 


s 


2 

Airplane wing area, ft 


t 


Time, seconds 


A t 


Increment of time, seconds 


V 

eg 


Airplane eg velocity, ft/sec 


X, Y.Z 


Airplane principal axis system 




Reference Cartesian axis system 


Z 

o 


Airplane altitude, feet 


Z 

o 


Rate of ascent, ft/sec 


oc 


Airplane angle of attack, alpha 


CU-. . 
IJ 


Direction cosines in terms of Euler angles (Table I) 


a 


Airplane sideslip angle, Beta 
Spin rate, per second 


^ a^<^ e 


Aileron, rudder and elevator deflections, degrees 


e.T.O 


Airplane principal axis system orientation Euler angles 
with respect to the 7 j ^ system 
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1 

f 




O Y ^ Aix'plane orientation in pitch, roll and yaw with respect 
^ to the Xj^, Yj, system 

f 



W 



Air density, slugs/ft' 
Aircraft weight, lbs 
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1 . INTRODUCTION 



A. BACKGROUND 

Aircraft flight tests frequently reveal significant differences between 
the predicted and the actual spin characteristics of high performance 
systems. Such differences often result in expensive engineering changes 
to the aircraft design and tend to reduce the mission capability of the 
aircraft. Furthermore, if the final spin characteristics are such that 
pilots are not permitted to practice the recovery from intentional 
spins, then many will consider the outer regions of the performance 
envelope with trepidation resulting in further degradation of the weapon 
system*s effectiveness. The need therefore, to accurately predict the 
spin characteristics of proposed aircraft designs is of paramount 
importance. 

Historically, efforts towards spin performance prediction have been 
concentrated in two areas: the development of system time histories 

and the analysis of free -spinning model test results. System time 
histories involve the simultaneous solution of the aircraft equations of 
motion with the required force and moment coefficients being deter- 
mined from wind tunnel model data. Solutions are obtained using 
digital computer programs which numerically integrate the appropriate 
equations of motion. There has been minimal success towards 
achieving a time history method that will accurately predict spin 
performance due to the critical dependence upon entry conditions. 
Free-spinning model tests utilize instrumented dynamically scaled 
models which are generally deployed into vertical wind tunnels. This 
method however has a limited capability in simulating the full scale 



11 



spin environment (Reynolds, Mach, and Fronde Numbers and the spin 
entry mechanism) and often the tests fail to identify all of the possible 
spin modes that may be experienced by the full scale vehicle. 

B, GENERAL 

This paper represents the final report on the development of a 
digital computer program designed to analytically predict the airplane 
spin .pe rformance using a unique set of aircraft equations of motion 
and simplifying assumptions regarding the nature of spin mode 
solutions. The primary hypothesis upon which the development of 
the computer program was predicated was the unique dynamically 
stable equilibrium (steady-state) solution for each aircraft spin mode. 
This assumption implies that all motion will ultimately be damped to 
a steady-state solution. Using the above hypothesis, the dynamically 
stable equilibrium solution may be determined directly without the 
encumbrances of a time history approach, with its required dependence 
upon entry conditions. The hypothesis ignores, however, the possible 
(and perhaps most probable) existence of a dynamically unstable 
equilibrium solution which could also yield a stable trajectory in 
*state spaceL 

C. SCOPE 

A rigorous research effort into the development of a means of 
predicting airplane spin modes would involve the following stages: 

I. Literature search 

II. Derivation of the mathematical model 

III. Verification of the mathematical model by using 
time history comparisons with baseline data 

IV. Development of a computer program designed to 
seek spin solutions 



12 



V. Verification of the computer program by utilizing 
wind tunnel data modified for simulator use 

VI. Determine and/or predict ’’real world” spin modes. 

Stages I and II are described in References [?] and [l]; while the 

base line data required for phase III is included in Reference [2]. The 
purpose of this effort was an attempt to accomplish steps IV and VI. 

This rather large undertaking was conducted concurrently with LCDR 
Buehler, who developed the applicable equations and coordinate systems 
required by stage II. The limited time available precluded a numerical 
validation of the computer program (phase III). 
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II. ANALYTIC CONSIDERATIONS 



A. COORDINATE SYSTEM 

The coordinate system used in the mathematical model is a com- 
bination of the airplane principal axis system and a fixed cylindrical 
system. The vertical (z) axis of the fixed cylindrical system repre- 
semts the initial central axis of motion. The cylindrical system 
locates the aircraft center of gravity in terms of the altitude coordinate 
(z^), the spin radius coordinate (R) and the angular position coordinate 
(C'). The orientation of the aircraft with respect to the cylindrical 
system is described in terms of Euler angles. A cartesian coordinate 
system is fixed at the eg position on the R vector with its axis (Xj ) in 
the (-R) direction and the axis is oriented in the (-z) direction. 

The cartesian system provides a reference for the orientation Euler 
angles whereby zero values of (0 ), and (0) would yield an upright, 

wings level airplane with its principal (x) axis pointed inward in the 
(-R) direction. Figure 1 depicts the details of the coordinate system; 
the positive directions of the position coordinates are also shown. 

As Buehler noted ’’The rationale behind this choice of coordinate 
system is that it more simply (in a mathematical sense) represents 
the motion being modeled. " [Reference l]. The advantage of this 
hybrid coordinate system is that it obviates a steady-state spin 
solution. The motion of the spinning aircraft is easily recognized by 
constant values of the six independent coordinate variables (R, z^, ^ , 

). Additionally the complexities of the stall, departure, and 
incipient spin motion have no impact on the results since the objective 
is to determine the steady-state spin modes. 
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FIXED REF, 



Fig. 1. The Coordinate System Used to Model the Spin 




Fig. 2 



A General Schematic of Euler Angle Relationships 



B. EULER ORIENTATION ANGLES 



The Euler Angle relationships utilized in the naathematical model 
involve rotations which are measured relative to the vertical (z^) axis 
and to the intersection of the (x, y) plane and the (x^^, ) plane. This 

system yields equations for the direction cosines as listed in Table I. 



Table I: Direction Cosines in Terms of Euler Angles 



Cosines of Antics between X,Y,Z and AT|, 


• 


X 


Y 


: z 




Cji = cos ^ cos 

— sin sin ^ cos $ 


C2i = — sin ^ COS 

— cos sin ^ cos 6 


= sin 8 sin 


Yi 


“ cos^sinvJ' 

-F sin^i cos\^ cos 8 


~ sin ^ sin ^ 

4- cos ^ cos ^ cos 6 


C32 = sin $ cos ^ 


Zi 


^13 — sin $ sin 


^23 = 8 cos 


fl33 = cos 8 



Inherent in this orientation system, however, is the difficulty in 
visualizing the final orientation without employing sketches or digital 
graphics equipment. This can be readily appreciated by referring to 
the general schematic of the Euler angle relationships depicted in 
Figure 2. In order to alleviate the visualization problem, Buehler, 
inReference [l], derived relationships between a set of orientation 
angles with a sequence of pitch (^p)> .and then yaw (C^y) 

and the Euler angles depicted in Figure 2. ^ The equations for 



Hereafter, and in the computer program, the orientation angles will 
be referred to as the ^ordered^ Euler angles, while the orientation 
angles depicted in Figure 2 will be referred to as the ’reference* 
Eule r angles. 



16 



converting the ordered Euler angles to the reference Euler angles are: 





//IaT^ 



( 










( 2 ) 




( 3 ) 



The equations for obtaining (®p)> 



from the reference 



Euler angles are: 



V ^ SMfSiL^) ( I Has j ' ’ 



0s> — 7 ^^^ 



.A rUScosi'^) ! _ 






Cos 



( 5 ) 
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The above equations have been numerically verified and validated 
for quadrants I and IV (jf 90^). This is sufficient for an aircraft 
oriented in an upright spin, however, for the case of an inverted spin 
where the aircraft will encounter quadrants II and III, the above 
relationships are invalid, 

C, THE EQUATIONS OF MOTION 

The spin is an uncontrolled large angle, six degree of freedom 
motion experienced by an aircraft in the stalled aerodynamic region. 
Since the spinning phenomenon is a complex motion which is influenced 
by a host of nonlinear variables, the equations which describe that 
motion are both nonlinear and coupled. Buehler, in Reference [l], 
derived the aircraft spin equations of motion utilized in the computer 
program. The equations and modifications thereto are described 
below. 

1. The Full Equations 

The equations of motion were cast in a more useable form by 
solving each equation for the corresponding acceleration parameter. 
This form avoids obscuration of the physical nature of the problem and 
places the numerical analyst in a better position to monitor, trouble- 
shoot and interpret iterative solution results. 
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2. The Modified Equations 



The full equations were further inodified by recognizing that 
only dynan^ically stable (or quasi- stable ) equilibrium solutions were 
sought. Thus angular orientation rates , CP ) were equated to zero, 

which in turn eliminated a considerable number of terms. The R 
parameter was retained because it effects the values of and 

to the extent that when ignored, it seriously impeded the iterative 
solution effort. 

3. The Short Equations 

The modified equations were further simplified by equating the 
R parameter to zero. This form is useful when seeking purely steady 
state solutions using the more mathematically elegant solution methods 
such as the gradient search method described in Section H. 

The full, modified, and short equations are listed in Appendix I. 



D. ANCILLIARY RELATIONSHIPS AND DEPENDENT VARIABLES 
In order to utilize conventional aerodynamic force and moment 
coefficients in the equations of motion, certain ancilliary relationships 
were needed. As derived by Buehler, and as depicted in Figure 3: 
Angle of Attack: 



CL - ^ -t- oLsi 'kc 



o4// /f -f- OL/2 



Sideslip Angle; 



(■7) 




7^ A' 



(oL-Zt R 3 oLz2- /■ oCzS Zo 

( oL„ /C ^ 



( 8 ) 
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Relative Wind Velocity: 



Roll Rate: 

^ « ( 'V* + ^^) 0^13 + O Cos Cj) 

Pitch 'Rate : 

O = (i- *f)o(z2 " © 5 '" ^ 



Yaw Rate : 

^ = {4 ^!k)’<s3 





Figure 3 

Definition of Positive c< and /? Angles 
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E. AERODYNAMIC FORCE AND MOMENT COEFFICIENTS 

The general form of .the aerodynamic force and moment coefficients 
employed in the equations of motion are as follows: 



^ mo ^ mo ^ m me) 



^ - X 



y^3 



( 13 ) 



0(^a) 'f ( ^0 /■ 



-f 



( 14 ) 



where represent the angular deflections of the ailerons, 

rubber and elevator respectively. 

All of the above coefficients on the right side of the equality sign 
are aerodynamic force and moment derivatives as defined in Table II. 
They must be experimentally determined and tabulated for a wide 
range of cxi and /? vlaues. f(p, q, or r) are functions of roll, pitch 
and yaw rates. 

The coefficients in Table II utilize c as a common characteristic 
length, thus allowing the use of a single length parameter input to the 
computer program. Since and j3 are zero for the steady spin, 
aerodynamic derivatives with respect to these two parameters have not 
been included in Table II. 
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SUBSCRIPT j 


i“x,y , 
or z 


1 


2 


3 


4 


5 


6 


7 








5Cfi 

m) 




^Cf^i 


dir 








^(■§) 




PCm/ 

5(ci) 









Table II 

2 

Definition of the Aerodynamic Derivatives 

The validity of any spin solution results is dependent upon how well 
the aerodynamic forces and moments have been modeled. One needs 
only to visualize the different dynamic environment that exists between 
the static model in a wind tunnel and the actual airplane experiencing 
three dimensional stalled flow in a spin. Herein lies the greatest 
weakness in the analytic spin solution effort. The static wind tunnel 
data utilized to determine the above 42 coefficients provides, at best, 
a very crude model of the actual conditions encountered by a spinning 
aircraft. However, for wanfof better data, the conventional data was 
utilized in order to exercise the program and to seek a first-order 
approximation to a spin solution. 

^Aerodynamic derivatives with respect to ^ and P have not been 
included sincee< and are zero for the steady spin. 
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F. SOLVING THE EQUATIONS OF MOTION 

A casual examination of the equations of motion found in Appendix 
II reveals that they have the following characteristics: 

1. ordinary second order differential eqrxations 

2. extensively coupled 

3. highly nonlinear 

4. force and moment terms are a function of three-dimensional 
flow over the geometry of the airplane. 

The above characteristics preclude direct analytical integration and 
thus one is forced to seek a numerical solution scheme which will pro- 
vide the greatest insight into the physical problem and also provide 
some degree of uniform convergence. Additionally, from numerical 
considerations the extreme length of the equations dictate the need for 
a method which maintains numerical significance and does not require 
excessive computational time. 

Two iterative solution methods were selected; the integration method 
which was utilized for steady and quasi- steady state solutions and the 
gradient method for the purely steady-state solution efforts. In both 
instances, the set of equations could be visualized as a '^black box” 
function transform with independent variables as inputs and the equation 
residuals as outputs. 

G. THE INTEGRATION SOLUTION METHOD 

The integration method can be visualized as the incremental 
tracking of the trajectory of a multiple degree of freedom spring-mass- 
damper system. The incremental changes are generated by ”psuedo” 
integration of the aerodynamic force and moment inbalances manifested 
in the equation residuals. This set of incremental changes constitutes 
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the "trajectory’’. The psuedo integration is accomplished by taking 
advantage of the form of the equations whereby the residuals are 
essentially the accelerations of the six spin degrees of freedom. 
Application of a common increment of time to elementary Euler 
integration formulae allows the computation of a d'^mamically consistent 
set of incremental changes in the independent variables and thus gener- 
ates new values for the subsequent iteration. Essentially, the method 
duplicates, in a numerical sense, what the aircraft experiences under 
actual flight conditions. 

One of the major advantages of casting the equations of motion in 
a cylindrical reference frame, is the ease with which tumbling motion 
can be separated from the motion associated with a steady state spin. 
Constant values of Z^, R, , 0 , , and d) indicate a steady state 

solution, while non-zero values of O , , or indicate tumbling 

motion. The computational scheme discards any angular momentum 
inconsistent with a steady-state spin and thus on each iteration any 
tumbling motion is suppressed. Such motion, if retained, could 
possibly carry the solution path (trajectory) through a steady state 
spin solution. 

The iterative procedure used in the computer was basically very 
simple. In response to a set of initial conditions (Zq, R, ^ , 

0 ) and related constants, the corresponding accelerations were 
computed. Application of a common time increment t ) to the basic 
Euler integration formulae yields a set of incremental changes which 
are then added to the previous set of independent variables to obtain a 
new set of initial conditions. The iteration process continues until the 
sum of the absolute values of the accelerations, denoted as S, is reduced 
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to less than some specified value. The ’’spring-mass-damper" system 
thus moves closer to equilibrium as it converges to a steady state 
solution. 

* * e € 

The acceleration Zq and R can be computed directly from the inde- 

pendent and dependent variables. The other accelerations ( 3^ , © , r , 

• • 

0 ) ) however, must be computed by a looping routine since the associ- 
ated equations are coupled to all of the angular accelerations. The 
looping ceases once a consistent set of accelerations is obtained. 

The incremental changes are obtained using the following equations: 



z 


(A/) 


<1 




( 15) 




II 


+(K At)h 




(16) 


Kohi) 


= Rm 


^ imd A-t + R 




(17) 




11 


^ A t) C 




(18) 




/ 


•^(e At) d 




(19) 




- 






(20) 
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•^(iJ At)c 




(21) 








(22) 




— 


->-(($ Ar)-f 




(23) 
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Figure 4 



A Schematic Representation of the 
Integration Solution Method 
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d>M + ( At $ ^‘)f 



(24) 



d). 






The constants (a, b, c, d, e, and f) are solution convergence factors 
which can be used to weigh the solution "rate” of a particular parameter. 
Also, one should note that the (G)i (U^), and ((p) expressions are used 
only in conjunction with the full equations. A schematic representation 
of the iterative process is depicted in Figure 4. 

H. THE GRADIENT SOLUTION METHOD 

The existance of a dynamically stable equilibrium solution can be 

verified utilizing a set of equations which exclude all time dependent 

terms. A gradient solution method was developed to search for such 

solutions using the previously discussed short equations. The method 

utilized is commonly referred to as the optimal method of steepest 
4 

descent. 

Conceptually, the solution space can be visualized as a series of 
peaks and valleys with solutions being the floor of any valley whose 
altitude is zero. The method is initiated at some 'geographic' position 
(initial conditions) where the "altitude" is measured (criterion function 
value). If the "altitude" were non-zero, the slope (gradient) of the 
"hill" would be measured and then the independent variables would be 
incrementally "marched" down the hill. The process continues until 
the independent variables reach the bottom of the "valley". 



Reference [?], 



page 2 75. 
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Mathematically, the criterion function is defined as the sum of the 
absolute values of the equation residuals. A solution is obtained when 



<o 



S = 2 ; )S;I =0 



(25) 



where S. are the equation residuals computed from the given independent 
variable vector, X. (i = 1, 6), 

The rapid reduction of S to zero is accomplished by iteratively 
choosing an increment such that <^5 approaches -S where 

SS- 



This is achieved by choosing6x;a s the product of the criterion gradient 
and a gain constant K; 



Sxy ~ K dX/ 



(27) 



However, by noting that 



<fx/ = ^ k dx,) = "S 



(28) 



or 



ss . ( g)‘ = -S 



(29) 



K can be expressed as 







( 30 ) 
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The subsequent value of the element of the X vector becomes 



X/. . - ix 



(AJrJ) 



Caj) 



OJ) 



'i 2 J 



( 31 ) 



where the subscript n refers to a sequential iteration number. 

is then computed as a function of if the condition 

1 



'/V 



(32) 



is satisfied, then the increment vector is reapplied, 



XaX/^/ ^ 



(33) 



A/ 



and this 1-D search process is continued in the direction of the original 
gradient until 



^ >9 



(34) 



whereupon new values of the gradients are computed. Since the process 
avoids the computation of the partial derivatives whenever the condition 
expressed in equation (32) is satisfied, significant amounts of computer 
computational time is saved and hence, the name ‘optimal’ gradient 
solution method is achieved. Figure 5 includes a flow chart of this 
technique. 
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Figure 5 



A Schematic Representation of the 
Gradient Solution Method 
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III. AIRCRAFT SPIN SOLUTION COMPUTER PROGRAM 



A, GENERAL 

This section presents briefly a general description of the digital 
computer program designed to solve the spin equations of motion for 
dynamically stable equilibrium solutions, A listing of the program as 
well as specific usage data and flow charts is presented in Appendix 
II. A program hierarchy is included as Figure 6 and an explanation of 
the program output is described by Figure 7. 

B, PROGRAM USAGE 

The program reads the tabular force /moment coefficients and 
related aircraft constants and then commences to iteratively seek a 
solution based upon the first initial condition data set that is read. 
Upon termination of a particular solution, the program restarts by 
reading the next initial condition data set, 

C, INPUT DATA 

The user must provide the following information to the program: 

1. Type aircraft (F4, Fill, etc. ) 

2. Tabular values of the force and moment coefficients/ 
stability derivatives (as described in Table II) and the 
corresponding values of angle of attack and sideslip 
angle 

3. Aircraft mass 

4. Atmospheric density 

5. Mean aerodynamic chord 

6. Gravitational constant 

7. Total wing area 
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Figure 6 

Airplane Spin Solution Program Hierarchy 
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8. Principal moments of inertia 

9. Control deflections 

10. Initial conditions of the independent variables 
The user also must stipulate the solution method, equation form, 
solution criteria, time increment, output format and solution con- 
vergence factors. 

The solution convergence factors available for the integration 
method, allow the user to accelerate or decelerate the solution ’’rate’' 
of any of the six independent variables. For example, by equating a 
particular factor to zero, that variable will remain constant, thus 
reducing the number of degrees of system freedom. 

D. OUTPUT DATA 

The printed output format is similar for the integration and gradient 
methods. The following is printed: 

1. Listing of all input data under appropriate headings 

2. Iterative results (integration method) 

a. iteration number (Ml) 

b. the number of looping iterations required to obtain 
a consistent set of angular accelerations (M2) 

c. the independent variables with their associated first 
and second derivatives 

d. relative wind velocity, angle of attack, and sideslip 
angle 

e. sum of the absolute values of the accelerations 

3. Iterative results (gradient method) 

a. ite ration numbe r (Ml) 

b. the iteration number corresponding to the last 
gradient computation (M2) 
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INTEGRATION METHOD ITERATIVE OUTPUT 
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c. the six independent variables with their corresponding 
residuals, incremental values, and partial derivatives 

d. relative wind velocity, angle of attack, and sideslip 
angle 

e. criterion value; the sum of the absolute values of the 
re siduals 

f. K; the solution gain factor 

4. Summary output 

Upon normal termination of the integration method, the iterative 
output is summarized. The listing includes the following data corres- 
ponding to every twenty-five iterations: 

1. iteration number 

2. independent variables 

3. dependent variables (o^ , /9 , ) 

Also, the above data corresponding to the smallest sum of the 
absolute values of the accelerations is printed. 

In the event of an abnormal termination, the program will print an 
appropriate message which informs the user as to the reason for the 
termination as well as the values of the offending variable(s). 

E. PROGRAM MODULE DESCRIPTION 

The program consists of nine modules; one main and eight sub- 
routines as described below: 

1. ' MAIN 

This routine is basically an input/output device. Tabulated 
data, fixed constants, and initial conditions are read into common 
memory, angular units are converted from degrees to radians; the 
optional output format is established and the first page is printed. 

The optional equation form and solution method, being determined. 
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a call is then made on the appropriate solution subroutine for sub- 
sequent computation, 

2. SOLVIT 

The integration solution method described in Section II is pro- 
grammed into this subroutine. The iterative results are printed and 
appropriate summary data is presented on the last page of the printout, 

3. GAWAIN 

This subroutine incorporate s the optimized -gradient method as 
described in Section II. Its operation is similar to SOLVIT, however, 
no summary data is printed. 

4. EONS 

Subroutine EQNS includes all of the independent and dependent 
equations. Calls are made primarily from either GAWAIN or SOLVIT. 
The current set of initial conditions is then used to compute the 
direction cosines, angle of attack and sideslip angle. Subroutine ANGLE 
is called to determine whether Alpha and Beta are within the appropri- 
ate definition limits. If so, the force and moment coefficients are then 
determined by interpolating the appropriate matrices using subroutines 
INTER, SPLINl, and SPLICO. 

The residuals are computed using either the short, modified or 
long form of the equations. A looping routine is utilized with the latter 
forms to obtain a consistent set of angular accelerations. Finally, 
upon computation of the sum of the absolute values of the residuals, 
execution is returned to the calling program. 

5. EULER 

The equations of motion require that angular units be in radians 
and that the reference Euler orientation angles be utilized. Subroutine 
EULER is used to convert computer variables to the output form of 
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angular units of degrees a.nd ordered orientation angles. The sub- 
routine is utilized by MAIN, GAWAIN and SOLVIT in order to allow 
output of the more usable ordered angles, 

6. ANGLE 

Subroutine ANGLE is utilized primarily by EQNS to check 
whether alpha or beta has exceeded the upper or lower limits of data 
definition. If the limits were exceeded, the iterative solution effort 
would be terminated and control returned to MAIN, ANGLE also reduces 
the absolute magnitude of angular displacements to less than 2TT radians. 

7. FIX 

Subroutine FIX is used in conjunction with IBM supplied ERRSET 
to facilitate termination of an iterative solution in the event of an under- 
flow, overflow or divide check while control is in EQNS. Control is 
returned to MAIN rather than effecting program termination as is the 
usual case. 

8. INTER 

Subroutine INTER is utilized by EQNS to interpolate the tabular 
force and moment coefficient data for specified values of alpha and 
beta. INTER calls SPLINl for the cubic interpolation of each matrix 
or vector. 

9. SPLINl and SPLICO 

Subroutines SPLINl and SPLICO were written by Mr. Rod Kure 
and are maintained by the NPS Computer Cente r ‘Staff. Their primary 
function is to provide cubic interpolation of a given matrix or vector. 
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V. AIRPLANE SPIN EQUATION VALIDATION 



The airplane spin equations of motion were validated by establish- 
ing initial conditions for a power-off glide configuration and then 
applying them to the equations using the spin solution program. 

A small computer program, TRIM, was prepared to enable the 
selection of a set of dynamically stable initial conditions. The TRIM 
program flowchart, the program listing and the sample output is 
included in Appendix D. 

The condition of a power-off glide is mathematically described by 
three expressions: 

(o< -o') ~ ( 34 ) 



V/ — ~k ^ V ^ dO'^ +• (fot-©)] (33) 












(36) 



where w, c< , © , ^ , C^, C^, and are as depicted in Figure 8. 

The values for and are determined by transposing from the 
principal to the stability axis systems as indicated below: 

Cl = ie) Co^ oc 4- (Cv ^e) Sjaj u (3V) 



Co Sy'^oi 



(38) 



Given constant values for o< , w, and S , one can solve equations 
(34)- -(36) for the quantities , © , and V, Subsequently the variables 
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Figure 8 

Schematic of the Aerodynamic Parameters for an Airplane in a Power-off Glide 



i 



2o ~ ^ ( oi-od 



(39) 



/? = \'^cos -0;> (40) 

can then be determined. 

The variables Z^, R, and are thus established independent 
of the spin equations of motion but based upon the same aerodynamic 
force and moment model utilized by the spin solution program. Appli- 
cation of these initial conditions to the spin solution program containing 
Buehler*s spin equations resulted in an immediate convergence whereby 
the sum of the absolute values of accelerations (equation residuals) was 
equal to 0. 0182. Table III provides a listing of the pertinent data 
utilized to achieve the equation validation. 



Table III 

Pertinent Data Relating to Equation Validation 





z 

o 


R 


Op 






TRIM Input: 


F-111 


force and moment model; air 


-2 

density = . 1553 x 10 


TRIM 

Output 


77. 01 


^ -381. d 


. 5763^ 


12. 0 - 


388.8 -5.534^ 


Integration 
Method 2 
Solution 


76. 99 


-381. 0 


. 57 


12. 0 


389. 0 



Notes: 1. Spin solution program input as per Table VIII 

2. SUM = . 0188 
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The spin solution program requires special care when exercising 
it for straight and level flight. Due to EUDER restrictions about 
accepting quadrant II and III and cardinal orientation angles, the 
initial condition values should be applied as in Table IV below. 

TABLE IV 

Initial Conditions for Equation Validation 



VARIABLE 


VALUE 




As used in TRIM Program 




As obtained from TRIM Output 




0. 0 




. 0001 


NS FLAG 


0 



The results of the equation validation computer tests conclusively 
validate the spin equations of motion for the power-off glide configu- 
ration. This procedure can only be considered, at best, qualitative, 
for it does not test the equations in a time dependent environment and 
configuration. Complete validation would require the generation of 
time histories and subsequent comparison to base line data such as is 
included in Reference [2]. 
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V. NUMERICAL ANALYSIS AND CONCLUSIONS 



A. AERODYNAMIC FORCE AND MOMENT MODELS 

The primary purpose of this thesis effort was to develop an aircraft 
spin solution computer program. Daring the program development, there 
were no idealized data available with known spin solutions, therefore the 
author utilized actual static wind tunnel data to provide an aerodynamic 
force and moment model. This was done in one case with the objective 
of ultimately comparing actual aircraft spin performance with computer 
generated predictions. 

Two models were utilized. Tabular F-111 data was obtained from 
Reference [3] and F-4b data were extracted from curves presented in 
Reference [4]. In both cases appropriate values were multiplied by 
b/c in order to normalize to the common dimension utilized in the 
computer program, i. e., the mean aerodynamic chord. The above 
data are included in Appendix E. 

B. SOLUTION METHOD EVALUATION 

1 . Integration Solution Method 

The integration solution method was the primary vehicle 
utilized to develop the spin equations of motion. As discussed in 
Section II, the primary advantage lies with the characteristic of the 
equation residuals being equivalent to the appropriate acceleration 
terms. 

The method must be initiated with a set of initial conditions. 

Many sets of various initial conditions were utilized during the 
development but for the sake of brevity only the results from one 
arbitrary set will be discussed. 
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It was evident early in the development phase that the modified 
equations provided the greatest potential. The full equations were 
suitable only for verifying the stability of solutions previously obtained 
from the modified equations. When the full equations were used to 
evaluate a set of arbitrary initial conditions, the angular orientation 
rates rapidly increased to disproportionate values and induced 
tumbling motion which, in turn, resulted in a divergent solution. The 
modified equations, however, enable quasi-steady solutions and are 
responsive to small changes in the time increment and control 
defle ctions . 

The solution paths generated by the integration method were 
characterized by an initially rapid convergence to a minimum value of 
the criterion function. The orientation values reorient to a nearly 
optimal configuration; the sink rate seeks the acceleration minimum; 
and the spin rate and radius effect trade-offs to achieve a consistent 
velocity value. Occasionally, the solution path became divergent in 
that the program was unable to find a consistent set of coupled 
accelerations. This was evidenced by the printout of the number of 
iterations required by the acceleration ”do loop” within Subroutine 
EQNS. The ”do loop” limit is set at twenty iterations to minimize 
computer time, however, only two iterations are usually needed for 
a stable solution path. Divergence is characterized by large and 
erratic changes in the independent variables such that alpha or beta 
exceed their respective fields of definition, thus resulting in solution 
termination. Solution path divergence can be minimized through use 
of a small time increment. The author determined that a range 
between . 001 and . 05 seconds proved optimal. The solution conver- 
gence factors were not utilized sufficiently to merit discussion. 
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Once the initial convergence was achieved, the individual 
solution paths seemed to vary in a seemingly arbitrary manner. 

While one variable tended toward a minimum acceleration value, other 
variables were producing increasing values. The variables continued 
to ’’wander”; changing magnitude and sign until either computer time 
was exhausted or the alpha or beta limits were exceeded. 

All solution efforts using F-4 data with pro-spin control 
'‘^deflections resulted in solution 'paths which exceeded the sideslip ^angle 
limits. Inclusion of rotary balance data to compensate for differential 
sideslip angles along the fuselage might inhibit such yawing action. 

Utilization of the integration method requires that the user 
learn to iterate the output to optimize the achievement of a final 
solution. This process is as described below: 

a. An arbitrary set of initial conditions is applied to the program 
to initiate a 500 iteration solution effort. This set of initial conditions 
should be such so as to approximate the anticipated spin solution mode. 

b. The initial condition set which yields the best convergence 

# 

is modified by substituting the value for corresponding to the 

#•# 

minimum in the 500 iteration output. An alternative is to substitute 
the set of three orientation angles corresponding to the minimum 
absolute sum of the orientation angle accelerations. Similarly, R and 
can be substituted. The resultant set of revised initial conditions 
is then used to initiate subsequent solution efforts of shorter {150 
iterations) duration. 

c. The process is continued until a best solution is obtained. 

d. The final solution is tested for stability by using the full 
equations with a small time increment. 
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¥ 1 ? 




Figure 9 and Table V illustrate a typical solution effort using an 
arbitrary set of initial conditions. 

The integration method results lead to a number of conclusions 

a. The F-111 aerodynamic force and moment model pre- 
cluded finding stable equilibrium solutions. This can be attributed to 
any of four reasons: (1) there were no truly steady state solutions; 

only oscillatory modes exist; (2) the "right" set of initial conditions 
was never employed and therefore a steady state solution path was not 
intercepted; (3) the steady state modes were lightly damped and there- 
fore, may be very slow in convergence; the program wasn’t exercised 
long enough to find the solution; and (4) the integration solution method 
itself precludes uniform convergence of all degrees of freedom. This 
would be due to self-excitation of the non-linear coupled acceleration 
te rms. 

b« The F-4 data do not model the actual aircraft since all 
initial conditions yielded solution paths with divergent yaw rates. The 
actual aircraft does not exhibit this trait and, in fact, demonstrates 
several steady and quasi-steady spin modes. 

c. The method is cumbersome to use because it requires 
that the user iteratively adjust inputed initial conditions in order to 
achieve a solution. 

d. The radius rate terms incorporated in the modified 
equations enables the determination of quasi-steady state solutions 
using the F-111 data. 

e. The orientation rate terms in the full equations were 
not compatible with the course integration techniques of this method. 
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Figure 9 

Schematic of a Typical Integration Solution Effort 
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Results of a Typical Integration Solution Effort 
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NOTES: 1. Run #7. utilized the gradient solution method. 

2. The inputted initial conditions cor re spoinding to Run #6 were also 
applied to the long equations using the integration method but the 
solution path was diverging when computation was terminated. 

3. F-111 aerodynamic model data was utilized with full pro-spin control deflections. 



f. The computer program needs to be modified to incorporate 
a subroutine which will compute the spin radius as well as the cylindrical 
inertial radius as in the present version. 

2. Gradient Solution Method 

The gradient solution method and its derivatives from optimal 
control theory have, perhaps, the greatest potential in the area of 
steady-stare spin prediction. The method programmed by the author 
can only be considered crude when compared to some of the mo;i’e 
sophisticated methods which utilize hybrid applications. 

The results obtained from a typical F-111 gradient solution 
effort provide the basis for evaluation of this method as well as yield 
insight into the F-111 solution space. Table VI includes a summary 
of these results. One should note that there was an initial convergence 
of the angular variables in the first fifty iterations. After iteration 
#75, the solution stabilized but the criterion value alternated between 
27. 0 and 102. 0; this was a result of the alternating residuals of the 
pitch and yaw equations. The sensitivity of the solution is extreme 
when one considers that a 0. 4% change in the absolute sum of the 
independent variables results in a 2, 740% reduction in the absolute 
sum of the acceleration. Specifically, the alternating changes in the 
Oand (f) residuals are the result of incremental changes in the © 
and (t) values of only 1 0 radians! 

The solution convergence becomes static due to the character- 

3 7 

istically large gradients (10 -10 ) and resultant small incremental 
step sizes. All gradient solutions attempted, resulted in no significant 
change in the values for and R parameters. 

Two conclusions can be made from the gradient solution method 

results; 
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Summary of a F-111 Gradient Solution Effort 
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NOTES: 1. Data is representative from iteration #25 to termination at #350. 

2. The © residual alternated between 0. 24 and + 61.2. 

3. The <t residual alternated between 0. 39 and + 14. 8. 



a. The method was unable to find a steady state solution; 
however, in one case, it was able to effect a reduction of the criterion 
function value by four orders of magnitude. 

b. The F-111 solution path is highly sensitive to small 
changes in the independent variable values. 

C. ANALYSIS OF THE F-111 AND F-4 STEADY-STATE SOLUTION 
SPACES 

The integration and gradient solution methods were both used to 

analyze the F-111 and F-4 steady state solution spaces. Initial 

conditions were generated by a 'Criterion Function Search Program'. 

This program evaluated criterion function values from a course grid 

of 186, 624 sets of initial conditions for each aircraft. 

One can visualize the solution space as a geographical relief with 

the 'solutions' occurring at the valley floors of zero elevation. The 

search program essentially measured the 'altitudes' (criterion function 

values) in this area and found only a very few points whose elevations 

4 

were less than 10 . Specifically, in the case of the F-4, only eleven 

sets of initial conditions generated criterion function values less than 
4 

10 and only one set of F-111 initial conditions met this criteria. The 
numerical results of the criterion function search program are included 
in Table VII. 

These generated initial conditions were evaluated using the gradient 
and integration solution methods. The results are extremely signifi- 
cant because they reveal that it is highly unlikely that a dynamically 
stable equilibrium solution exists. These solution results are 
summarized in Figure 10 and Table VIII. 

The two conclusions which can be made are: 
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!• The static wind tunnel data utilized to form the basis for the 
F-111 and F-4 solution spaces are of such a nature as to preclude 
finding a dynamically stable equilibrium solution. There was no 
solution convergence to a criterion value of less than Z5. 0. 

2. The aerodynamic data utilized do not adequately describe the 
actual aerodynamic force and moment field experienced by the actual 
aircraft. This is evidenced by the failure to find even a quasi-steady 
solution similar to flight test results. 

The above implications are not surprising considering that the 
aerodynamics of the spin problem involve low speed, three-dimensional 
stalled flow and, as such, forces and moments are nonlinear. Thus 
the superposition of the separate effects of triaxial rotation, control 
deflection, and other static wind tunnel data cannot be completely 
justified analytically. The only wind tunnel test technique that 
presently appears to offer any possibility of yielding satisfactory 
aerodynamic data for the spin is the rotary balance force and moment 
measuring apparatus which is described in Reference [6]. 
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VI. SUMMARY OF THE MAJOR CONCH US IQ NS 



The criterion function search revealed that, for the F-111 and F-4 
aerodynamic force and moment models utilized in this study, it was 
highly improbable that a dynamically stable equilibrium solution 
existed. 

The airplane spin equations described in Reference [l] were 
qualitatively validated. Quantitative verification is needed and it is 
recommended that the equations be cast in a form compatible with a 
time history analysis in order that the resultant data may be compared 
with the base-line spin data included in Reference [2]. 

The integration solution method is too cumbersome to be effectively 
utilized although it can provide information relating to quasi-steady 
solutions once the user has learned to iteratively manipulate the initial 
conditions . 

The gradient solution method demonstrated its potential by effecting 
significant reductions in criterion function values to quasi- steady state 
solutions. 
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APPENDIX A 



FUEL SPIN EQUATIONS 

The following equations constit\ite the full eq^iations of motion for 
a spinning airplane as written in a cylindrical coordinate system. 



Z Eqriation 

\ 

/V ^ ^ (lfyO<23 ^"< 5 ^ 



R Equation 
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MODIFIED EQUATIONS 
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© Equation: 












Equation: 
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SHORT EQUATIONS 



Equation: 
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APPENDIX B 



AIRCRAFT SPIN SOLUTION PROGRAM 

A. GENERAL 

This appendix is intended for the Fortran programmer who is a 
potential user of the program. The included source listing contains 
an. extensive introductory comment section with which the reader 
should familiarize himself prior to reading the other sections of this 
appendix. 

B. DATA INPUT 

The composition of the data deck is described in the source listing 
and also depicted in Figure B-1. One should particularly note the 
composition of data sets ^2--#43, The first card of each of the afore- 
mentioned sets ’’tells^’ the program which tabular data is to follow, 
its literal name, and whether the data set is a matrix, a vector, or a 
null set. The subsequent cards of each data set contain the tabular data 
except where no data is available (null set); in which case, the data set 
consists only of the first card. Each force /moment coefficient matrix 
is characterized by a unique integer assignment (NA) as specified in 
Table B-I. NA is the first value of each of the 42 coefficient data sets. 
There must be exactly 42 ’’NA” cards however there is no restriction 
as to their sequence. 

Each solution effort is initiated based upon the initial conditions 
and associated program variables provided by data set #45. An option 
is provided (integration method only) for the program variables to 
remain constant for subsequent solution efforts; only the initial 
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Data Deck Composition for the Airplane Spin Solution Program 
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conditions included as data set #48 need be provided* This option is 
effected by assigning NSFLAG = 2 or 3 as explained in the program 
comments listing and Table B~II. 

It is suggested that the user who wishes to test a large number of 
initial conditions assign LIMIT = 10, NSFLAG = 2 and then test all of 
the initial conditions on one job. This scheme will quickly reveal to the 
user which initial conditions exceed the alpha/beta limits without 
wasting valuable computation time. Subsequently, those initial con- 
ditions which appear to have a solution potential can be submitted with 
appropriate data set #45 values. 

Since every iteration generates a new set of initial conditions, the 
program can be restarted merely by including as data set #45 the appropri- 
ate values taken from the printed output. This option allows the user to 
check the stability of a quasi-steady spin solution by changing the time 
increment (BELT), one of the solution convergence factors (SCF), the 
control deflections (AILDEF, RUDDEF, ELEDEF), or the atmospheric 
density (RHO). 

C. DATA OUTPUT 

All output is formated for 130 character computer paper. An option 
is provided by NPRINT to suppress the iterative results; otherwise, all 
input data is echo checked, all significant iterative computations are 
printed and a summary of the integration solution path is provided. 

D. STORAGE REQUIREMENTS 

The amount of storage required is primarily a function of the size 
of the CFORCM three dimensional matrix. The third index of the 
matrix is associated with, and is no larger than, the largest NA 
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TABLE B-II 



NSFLAG Integer Assignment 



NSFLAG 


SOLUTION 

METHOD 


EQUATION 

FORM 


SUBSEQUENT 
DATA SET 
COMPOSITION 


0 


Integration 


Modified 


#45 
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Integration 


Full 


#45 
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Integration 


Modified 


#48 


3 


Integration 


. Full 


#48 


4 


Gradient 


Short 


#45 
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integer (NSMAX) corresponding to vector or matrix coefficient data sets. 
The first and second indices are associated with IMAX and JMAX. 

Prior to using the program the user should dimension A, B, and 
CFORCM in the COMMON/ WORKA/ specification statements to 
A(IMAX), B(JMAX), and CFORCM(IMAX, JMAX, NAMAX). The 
WORKA statement must be identical in the following routines: MAIN, 

GAWAIN, SOLVIT, EQNS, INTER, and ANGLE, Common regions are 
utilized for most communications between subroutines and therefore 
care must be exercised to ensure that the common statements are 
identical. 

E. PROGRAM FLOWCHARTS 

Main and subroutine flow charts are included as Figures B-2 to B-9 
and are intended to be self-explanatory. The program source listing 
includes a sufficient number of comments to amplify the flowchart 
information. Subroutines SPLINl and SPLICO were provided by the 
NPS Computer Center. No effort has been made to provide reference 
flowcharts as both routines have been extensively tested and validated 
by the NPS staff. 

F. PROGRAMMING CONSIDERATIONS 

One should note that calls on Euler are made in pairs. The angular 
variables are stored in the common memory, thus if the first conver- 
sion is from equation compatible to output variables, the second call 
must effect the opposite conversion. 

Subroutine GAWAIN requires the computation of the partial derivative 
of the criterion with respect to a change in each independent variable. 
This is accomplished by incrementing a single independent variable by 
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FAC. The criterion is then computed using this value. The change in 
the criterion function value is then divided by FAC to obtain the partial 
de rivative . 
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COMPUTE ALPHA AND BETA 
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SOLVIT Subroutine Flowchart 
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GAWAIN Subroutine Flowchart 
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EQNS Subroutine Flowchart 



EULER (N) 



r COMMON VARIABLES UTILIZED ! 




Figure B-6 

EULER Subroutine Flowchart 
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S UBROU T IxNE ANGLE ( *) j 



COMMON VARIABLES UTILIZED i 

I 

A, B, ALPHA, BETA, PI, ! 

ALF, BET, R, NFLAG, IMAX, JAMX ^ ; 




Figure B-7 

ANGLE Subroutine Flowchart 
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f SUBROUTINE INTER (NA, VALUE) 




COMMON VARIABLES UTILIZED 



CFORCM NA KEY 

i A IMAX ALPHA 

I B JMAX __BETA_ 





Creturn^ 



Figure B-8 

INTER Subroutine Flowchart 
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SUBROUTINE . 


FIX (IRETCD) I 




1 
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1 
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1 

1 



1 COMMON VARIABLE: NFLAG 1 



NFLAG = 1 
IRETCD = 0 



(RETURN ~y 



Figure B-9 

FIX Subroutine Flowchart 
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APPENDIX C 



CRITERION FUNCTION SEARCH PROGRAM 
The aircraft spin solution program was modified in order to evaluate 
the "reasonableness'' of various sets of initial conditions. The result- 
ing Criterion Function Search program is designed to compute the 
criterion value for a given set of initial conditions. The initial con- 
ditions are generated by six nested "do loops" within MAIN. The 
dependent variables aC , /3 , V and n are subsequently computed and 

compared to pre-established constraints; if the constraints are satisfied, 

4 

the criterion is subsequently computed. If the value is less than 10 , the 
independent variables are printed. An abbreviated flowchart of the pro- 
gram is included as Figure C-1. Pertinent data relating to the initial 
condition "do loops" and dependent variable constraints is included as 
Table C-I. The program is basically the same as the spin solution 
program except that SOLVIT and GAWAIN are not included and that 
MAIN and EQNS are modified to include the initial condition "do loops" 
and constraints. The range of the do loops is such that the initial 
conditions are all within a reasonable approximation to an erect spin. 
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Figure C-1 

Abbreviated Flowcha.rt for the Criterion Function 
Search Computer Program 
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TABLE C-I 



Pertinent Data Regarding the Criterion Function 
Search Program Independent and Dependent Variables 



Initial Condition *'Do Loop^' Ranges 



c 

Z. 


R 




©p 




CD 


-150. 0 


25. 0 


10. 0 


-11. 0 


5. 0 


- 5. 


-200. 0 


100. 0 


20. 0 


-21. 0 


15, 0 


-15, 


-250. 0 


175. 0 


30, 0 


-31. 0 


25. 0 


-25. 


-300. 0 


250. 0 


40. 0 


-41. 0 


35. 0 


-35. 


-350. 0 


325. 0 


50. 0 


-51. 0 


45. 0 


-45. 


-400. 0 


400. 0 


60. 0 


-61, 0 


55, 0 


-55, 

-65. 

-75. 

-85. 



Dependent Variable Constraints 

= 0° + 40° 

= -1. 0° to + 89. 0° 

^ 589. 0 (corresponds to 300 knots equivalent airspeed) 
^ 2.0 g's 



Fixed Constants 



II 


. 001756 


slugs 






F-111 




F-4 


= 


10. 0 




21. 


II 


-15. 0 




- 30, 


Sg = 


-25. 0 




- 30. 



y 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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APPENDIX E 



AERODYNAMIC FORCE AND MOMENT MODELS 



F-111 Data 


A. Mass 


1555. 0 


slugs 


Chord 


9. 04 


feet 


Wing Area 


525 


feet^ 


I 

X 


50, 000 


slug ft' 


I 


315, 200 


slug ft' 


y 


I 

Z 


351, 500 


slug ft' 


B. Attached data F- 


111 reduced to 6 x 


F-4 Data 


A. Mass 


1345 


slugs 


Chord 


16. 03 


feet 


Wing Area 


530 


feet^ 


I 

X 


23, 000 


slug ft 


I 


134, 000 


slug ft 


y 


I 

z 


152, 000 


slug ft 


B. Attached data (F 


-4) reduced 


to 6 X 
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